1//////////////////////////////////////////////////////////////////////
2// LibFile: math.scad
3// Assorted math functions, including linear interpolation, list operations (sums, mean, products),
4// convolution, quantization, log2, hyperbolic trig functions, random numbers, derivatives,
5// polynomials, and root finding.
6// Includes:
7// include <BOSL2/std.scad>
8// FileGroup: Math
9// FileSummary: Math on lists, special functions, quantization, random numbers, calculus, root finding
10//
11// FileFootnotes: STD=Included in std.scad
12//////////////////////////////////////////////////////////////////////
13
14// Section: Math Constants
15
16// Constant: PHI
17// Synopsis: The golden ratio φ (phi). Approximately 1.6180339887
18// Topics: Constants, Math
19// See Also: EPSILON, INF, NAN
20// Description: The golden ratio φ (phi). Approximately 1.6180339887
21PHI = (1+sqrt(5))/2;
22
23// Constant: EPSILON
24// Synopsis: A tiny value to compare floating point values. `1e-9`
25// Topics: Constants, Math
26// See Also: PHI, EPSILON, INF, NAN
27// Description: A really small value useful in comparing floating point numbers. ie: abs(a-b)<EPSILON `1e-9`
28EPSILON = 1e-9;
29
30// Constant: INF
31// Synopsis: The floating point value for Infinite.
32// Topics: Constants, Math
33// See Also: PHI, EPSILON, INF, NAN
34// Description: The value `inf`, useful for comparisons.
35INF = 1/0;
36
37// Constant: NAN
38// Synopsis: The floating point value for Not a Number.
39// Topics: Constants, Math
40// See Also: PHI, EPSILON, INF, NAN
41// Description: The value `nan`, useful for comparisons.
42NAN = acos(2);
43
44
45
46// Section: Interpolation and Counting
47
48
49// Function: count()
50// Synopsis: Creates a list of incrementing numbers.
51// Topics: Math, Indexing
52// See Also: idx()
53// Usage:
54// list = count(n, [s], [step], [reverse]);
55// Description:
56// Creates a list of `n` numbers, starting at `s`, incrementing by `step` each time.
57// You can also pass a list for n and then the length of the input list is used.
58// Arguments:
59// n = The length of the list of numbers to create, or a list to match the length of
60// s = The starting value of the list of numbers.
61// step = The amount to increment successive numbers in the list.
62// reverse = Reverse the list. Default: false.
63// Example:
64// nl1 = count(5); // Returns: [0,1,2,3,4]
65// nl2 = count(5,3); // Returns: [3,4,5,6,7]
66// nl3 = count(4,3,2); // Returns: [3,5,7,9]
67// nl4 = count(5,reverse=true); // Returns: [4,3,2,1,0]
68// nl5 = count(5,3,reverse=true); // Returns: [7,6,5,4,3]
69function count(n,s=0,step=1,reverse=false) = let(n=is_list(n) ? len(n) : n)
70 reverse? [for (i=[n-1:-1:0]) s+i*step]
71 : [for (i=[0:1:n-1]) s+i*step];
72
73
74// Function: lerp()
75// Synopsis: Linearly interpolates between two values.
76// Topics: Interpolation, Math
77// See Also: v_lookup(), lerpn()
78// Usage:
79// x = lerp(a, b, u);
80// l = lerp(a, b, LIST);
81// Description:
82// Interpolate between two values or vectors.
83// If `u` is given as a number, returns the single interpolated value.
84// If `u` is 0.0, then the value of `a` is returned.
85// If `u` is 1.0, then the value of `b` is returned.
86// If `u` is a range, or list of numbers, returns a list of interpolated values.
87// It is valid to use a `u` value outside the range 0 to 1. The result will be an extrapolation
88// along the slope formed by `a` and `b`.
89// Arguments:
90// a = First value or vector.
91// b = Second value or vector.
92// u = The proportion from `a` to `b` to calculate. Standard range is 0.0 to 1.0, inclusive. If given as a list or range of values, returns a list of results.
93// Example:
94// x = lerp(0,20,0.3); // Returns: 6
95// x = lerp(0,20,0.8); // Returns: 16
96// x = lerp(0,20,-0.1); // Returns: -2
97// x = lerp(0,20,1.1); // Returns: 22
98// p = lerp([0,0],[20,10],0.25); // Returns [5,2.5]
99// l = lerp(0,20,[0.4,0.6]); // Returns: [8,12]
100// l = lerp(0,20,[0.25:0.25:0.75]); // Returns: [5,10,15]
101// Example(2D):
102// p1 = [-50,-20]; p2 = [50,30];
103// stroke([p1,p2]);
104// pts = lerp(p1, p2, [0:1/8:1]);
105// // Points colored in ROYGBIV order.
106// rainbow(pts) translate($item) circle(d=3,$fn=8);
107function lerp(a,b,u) =
108 assert(same_shape(a,b), "Bad or inconsistent inputs to lerp")
109 is_finite(u)? (1-u)*a + u*b :
110 assert(is_finite(u) || is_vector(u) || valid_range(u), "Input u to lerp must be a number, vector, or valid range.")
111 [for (v = u) (1-v)*a + v*b ];
112
113
114// Function: lerpn()
115// Synopsis: Returns exactly `n` values, linearly interpolated between `a` and `b`.
116// Topics: Interpolation, Math
117// See Also: v_lookup(), lerp()
118// Usage:
119// x = lerpn(a, b, n);
120// x = lerpn(a, b, n, [endpoint]);
121// Description:
122// Returns exactly `n` values, linearly interpolated between `a` and `b`.
123// If `endpoint` is true, then the last value will exactly equal `b`.
124// If `endpoint` is false, then the last value will be `a+(b-a)*(1-1/n)`.
125// Arguments:
126// a = First value or vector.
127// b = Second value or vector.
128// n = The number of values to return.
129// endpoint = If true, the last value will be exactly `b`. If false, the last value will be one step less.
130// Example:
131// l = lerpn(-4,4,9); // Returns: [-4,-3,-2,-1,0,1,2,3,4]
132// l = lerpn(-4,4,8,false); // Returns: [-4,-3,-2,-1,0,1,2,3]
133// l = lerpn(0,1,6); // Returns: [0, 0.2, 0.4, 0.6, 0.8, 1]
134// l = lerpn(0,1,5,false); // Returns: [0, 0.2, 0.4, 0.6, 0.8]
135function lerpn(a,b,n,endpoint=true) =
136 assert(same_shape(a,b), "Bad or inconsistent inputs to lerpn")
137 assert(is_int(n))
138 assert(is_bool(endpoint))
139 let( d = n - (endpoint? 1 : 0) )
140 [for (i=[0:1:n-1]) let(u=i/d) (1-u)*a + u*b];
141
142
143
144// Section: Miscellaneous Functions
145
146// Function: sqr()
147// Synopsis: Returns the square of the given value.
148// Topics: Math
149// See Also: hypot(), log2()
150// Usage:
151// x2 = sqr(x);
152// Description:
153// If given a number, returns the square of that number,
154// If given a vector, returns the sum-of-squares/dot product of the vector elements.
155// If given a matrix, returns the matrix multiplication of the matrix with itself.
156// Example:
157// sqr(3); // Returns: 9
158// sqr(-4); // Returns: 16
159// sqr([2,3,4]); // Returns: 29
160// sqr([[1,2],[3,4]]); // Returns [[7,10],[15,22]]
161function sqr(x) =
162 assert(is_finite(x) || is_vector(x) || is_matrix(x), "Input is not a number nor a list of numbers.")
163 x*x;
164
165
166// Function: log2()
167// Synopsis: Returns the log base 2 of the given value.
168// Topics: Math
169// See Also: hypot(), sqr()
170// Usage:
171// val = log2(x);
172// Description:
173// Returns the logarithm base 2 of the value given.
174// Example:
175// log2(0.125); // Returns: -3
176// log2(16); // Returns: 4
177// log2(256); // Returns: 8
178function log2(x) =
179 assert( is_finite(x), "Input is not a number.")
180 ln(x)/ln(2);
181
182// this may return NAN or INF; should it check x>0 ?
183
184// Function: hypot()
185// Synopsis: Returns the hypotenuse length of a 2D or 3D triangle.
186// Topics: Math
187// See Also: hypot(), sqr(), log2()
188// Usage:
189// l = hypot(x, y, [z]);
190// Description:
191// Calculate hypotenuse length of a 2D or 3D triangle.
192// Arguments:
193// x = Length on the X axis.
194// y = Length on the Y axis.
195// z = Length on the Z axis. Optional.
196// Example:
197// l = hypot(3,4); // Returns: 5
198// l = hypot(3,4,5); // Returns: ~7.0710678119
199function hypot(x,y,z=0) =
200 assert( is_vector([x,y,z]), "Improper number(s).")
201 norm([x,y,z]);
202
203
204// Function: factorial()
205// Synopsis: Returns the factorial of the given integer.
206// Topics: Math
207// See Also: hypot(), sqr(), log2()
208// Usage:
209// x = factorial(n, [d]);
210// Description:
211// Returns the factorial of the given integer value, or n!/d! if d is given.
212// Arguments:
213// n = The integer number to get the factorial of. (n!)
214// d = If given, the returned value will be (n! / d!)
215// Example:
216// x = factorial(4); // Returns: 24
217// y = factorial(6); // Returns: 720
218// z = factorial(9); // Returns: 362880
219function factorial(n,d=0) =
220 assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is defined only for non negative integers")
221 assert(d<=n, "d cannot be larger than n")
222 product([1,for (i=[n:-1:d+1]) i]);
223
224
225// Function: binomial()
226// Synopsis: Returns the binomial coefficients of the integer `n`.
227// Topics: Math
228// See Also: hypot(), sqr(), log2(), factorial()
229// Usage:
230// x = binomial(n);
231// Description:
232// Returns the binomial coefficients of the integer `n`.
233// Arguments:
234// n = The integer to get the binomial coefficients of
235// Example:
236// x = binomial(3); // Returns: [1,3,3,1]
237// y = binomial(4); // Returns: [1,4,6,4,1]
238// z = binomial(6); // Returns: [1,6,15,20,15,6,1]
239function binomial(n) =
240 assert( is_int(n) && n>0, "Input is not an integer greater than 0.")
241 [for( c = 1, i = 0;
242 i<=n;
243 c = c*(n-i)/(i+1), i = i+1
244 ) c ] ;
245
246
247// Function: binomial_coefficient()
248// Synopsis: Returns the `k`-th binomial coefficient of the integer `n`.
249// Topics: Math
250// See Also: hypot(), sqr(), log2(), factorial()
251// Usage:
252// x = binomial_coefficient(n, k);
253// Description:
254// Returns the `k`-th binomial coefficient of the integer `n`.
255// Arguments:
256// n = The integer to get the binomial coefficient of
257// k = The binomial coefficient index
258// Example:
259// x = binomial_coefficient(3,2); // Returns: 3
260// y = binomial_coefficient(10,6); // Returns: 210
261function binomial_coefficient(n,k) =
262 assert( is_int(n) && is_int(k), "Some input is not a number.")
263 k < 0 || k > n ? 0 :
264 k ==0 || k ==n ? 1 :
265 let( k = min(k, n-k),
266 b = [for( c = 1, i = 0;
267 i<=k;
268 c = c*(n-i)/(i+1), i = i+1
269 ) c] )
270 b[len(b)-1];
271
272
273// Function: gcd()
274// Synopsis: Returns the Greatest Common Divisor/Factor of two integers.
275// Topics: Math
276// See Also: hypot(), sqr(), log2(), factorial(), binomial(), gcd(), lcm()
277// Usage:
278// x = gcd(a,b)
279// Description:
280// Computes the Greatest Common Divisor/Factor of `a` and `b`.
281function gcd(a,b) =
282 assert(is_int(a) && is_int(b),"Arguments to gcd must be integers")
283 b==0 ? abs(a) : gcd(b,a % b);
284
285
286// Computes lcm for two integers
287function _lcm(a,b) =
288 assert(is_int(a) && is_int(b), "Invalid non-integer parameters to lcm")
289 assert(a!=0 && b!=0, "Arguments to lcm should not be zero")
290 abs(a*b) / gcd(a,b);
291
292
293// Computes lcm for a list of values
294function _lcmlist(a) =
295 len(a)==1 ? a[0] :
296 _lcmlist(concat(lcm(a[0],a[1]),list_tail(a,2)));
297
298
299// Function: lcm()
300// Synopsis: Returns the Least Common Multiple of two or more integers.
301// Topics: Math
302// See Also: hypot(), sqr(), log2(), factorial(), binomial(), gcd(), lcm()
303// Usage:
304// div = lcm(a, b);
305// divs = lcm(list);
306// Description:
307// Computes the Least Common Multiple of the two arguments or a list of arguments. Inputs should
308// be non-zero integers. The output is always a positive integer. It is an error to pass zero
309// as an argument.
310function lcm(a,b=[]) =
311 !is_list(a) && !is_list(b)
312 ? _lcm(a,b)
313 : let( arglist = concat(force_list(a),force_list(b)) )
314 assert(len(arglist)>0, "Invalid call to lcm with empty list(s)")
315 _lcmlist(arglist);
316
317// Function rational_approx()
318// Usage:
319// pq = rational_approx(x, maxq);
320// Description:
321// Finds the best rational approximation p/q to the number x so that q<=maxq. Returns
322// the result as `[p,q]`. If the input is zero, then returns `[0,1]`.
323// Example:
324// pq1 = rational_approx(PI,10); // Returns: [22,7]
325// pq2 = rational_approx(PI,10000); // Returns: [355, 113]
326// pq3 = rational_approx(221/323,500); // Returns: [13,19]
327// pq4 = rational_approx(0,50); // Returns: [0,1]
328function rational_approx(x, maxq, cfrac=[], p, q) =
329 let(
330 next = floor(x),
331 fracpart = x-next,
332 cfrac = [each cfrac, next],
333 pq = _cfrac_to_pq(cfrac)
334 )
335 approx(fracpart,0) ? pq
336 : pq[1]>maxq ? [p,q]
337 : rational_approx(1/fracpart,maxq,cfrac, pq[0], pq[1]);
338
339
340// Converts a continued fraction given as a list with leading integer term
341// into a fraction in the form p / q, returning [p,q].
342function _cfrac_to_pq(cfrac,p=0,q=1,ind) =
343 is_undef(ind) ? _cfrac_to_pq(cfrac,p,q,len(cfrac)-1)
344 : ind==0 ? [p+q*cfrac[0], q]
345 : _cfrac_to_pq(cfrac, q, cfrac[ind]*q+p, ind-1);
346
347
348// Section: Hyperbolic Trigonometry
349
350// Function: sinh()
351// Synopsis: Returns the hyperbolic sine of the given value.
352// Topics: Math, Trigonometry
353// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
354// Usage:
355// a = sinh(x);
356// Description: Takes a value `x`, and returns the hyperbolic sine of it.
357function sinh(x) =
358 assert(is_finite(x), "The input must be a finite number.")
359 (exp(x)-exp(-x))/2;
360
361
362// Function: cosh()
363// Synopsis: Returns the hyperbolic cosine of the given value.
364// Topics: Math, Trigonometry
365// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
366// Usage:
367// a = cosh(x);
368// Description: Takes a value `x`, and returns the hyperbolic cosine of it.
369function cosh(x) =
370 assert(is_finite(x), "The input must be a finite number.")
371 (exp(x)+exp(-x))/2;
372
373
374// Function: tanh()
375// Synopsis: Returns the hyperbolic tangent of the given value.
376// Topics: Math, Trigonometry
377// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
378// Usage:
379// a = tanh(x);
380// Description: Takes a value `x`, and returns the hyperbolic tangent of it.
381function tanh(x) =
382 assert(is_finite(x), "The input must be a finite number.")
383 sinh(x)/cosh(x);
384
385
386// Function: asinh()
387// Synopsis: Returns the hyperbolic arc-sine of the given value.
388// Topics: Math, Trigonometry
389// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
390// Usage:
391// a = asinh(x);
392// Description: Takes a value `x`, and returns the inverse hyperbolic sine of it.
393function asinh(x) =
394 assert(is_finite(x), "The input must be a finite number.")
395 ln(x+sqrt(x*x+1));
396
397
398// Function: acosh()
399// Synopsis: Returns the hyperbolic arc-cosine of the given value.
400// Topics: Math, Trigonometry
401// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
402// Usage:
403// a = acosh(x);
404// Description: Takes a value `x`, and returns the inverse hyperbolic cosine of it.
405function acosh(x) =
406 assert(is_finite(x), "The input must be a finite number.")
407 ln(x+sqrt(x*x-1));
408
409
410// Function: atanh()
411// Synopsis: Returns the hyperbolic arc-tangent of the given value.
412// Topics: Math, Trigonometry
413// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
414// Usage:
415// a = atanh(x);
416// Description: Takes a value `x`, and returns the inverse hyperbolic tangent of it.
417function atanh(x) =
418 assert(is_finite(x), "The input must be a finite number.")
419 ln((1+x)/(1-x))/2;
420
421
422// Section: Quantization
423
424// Function: quant()
425// Synopsis: Returns `x` quantized to the nearest integer multiple of `y`.
426// Topics: Math, Quantization
427// See Also: quant(), quantdn(), quantup()
428// Usage:
429// num = quant(x, y);
430// Description:
431// Quantize a value `x` to an integer multiple of `y`, rounding to the nearest multiple.
432// The value of `y` does NOT have to be an integer. If `x` is a list, then every item
433// in that list will be recursively quantized.
434// Arguments:
435// x = The value or list to quantize.
436// y = Positive quantum to quantize to
437// Example:
438// a = quant(12,4); // Returns: 12
439// b = quant(13,4); // Returns: 12
440// c = quant(13.1,4); // Returns: 12
441// d = quant(14,4); // Returns: 16
442// e = quant(14.1,4); // Returns: 16
443// f = quant(15,4); // Returns: 16
444// g = quant(16,4); // Returns: 16
445// h = quant(9,3); // Returns: 9
446// i = quant(10,3); // Returns: 9
447// j = quant(10.4,3); // Returns: 9
448// k = quant(10.5,3); // Returns: 12
449// l = quant(11,3); // Returns: 12
450// m = quant(12,3); // Returns: 12
451// n = quant(11,2.5); // Returns: 10
452// o = quant(12,2.5); // Returns: 12.5
453// p = quant([12,13,13.1,14,14.1,15,16],4); // Returns: [12,12,12,16,16,16,16]
454// q = quant([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,12,12,12]
455// r = quant([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[12,12,12]]
456function quant(x,y) =
457 assert( is_finite(y) && y>0, "The quantum `y` must be a positive value.")
458 is_num(x) ? round(x/y)*y
459 : _roundall(x/y)*y;
460
461function _roundall(data) =
462 [for(x=data) is_list(x) ? _roundall(x) : round(x)];
463
464
465// Function: quantdn()
466// Synopsis: Returns `x` quantized down to an integer multiple of `y`.
467// Topics: Math, Quantization
468// See Also: quant(), quantdn(), quantup()
469// Usage:
470// num = quantdn(x, y);
471// Description:
472// Quantize a value `x` to an integer multiple of `y`, rounding down to the previous multiple.
473// The value of `y` does NOT have to be an integer. If `x` is a list, then every item in that
474// list will be recursively quantized down.
475// Arguments:
476// x = The value or list to quantize.
477// y = Postive quantum to quantize to.
478// Example:
479// a = quantdn(12,4); // Returns: 12
480// b = quantdn(13,4); // Returns: 12
481// c = quantdn(13.1,4); // Returns: 12
482// d = quantdn(14,4); // Returns: 12
483// e = quantdn(14.1,4); // Returns: 12
484// f = quantdn(15,4); // Returns: 12
485// g = quantdn(16,4); // Returns: 16
486// h = quantdn(9,3); // Returns: 9
487// i = quantdn(10,3); // Returns: 9
488// j = quantdn(10.4,3); // Returns: 9
489// k = quantdn(10.5,3); // Returns: 9
490// l = quantdn(11,3); // Returns: 9
491// m = quantdn(12,3); // Returns: 12
492// n = quantdn(11,2.5); // Returns: 10
493// o = quantdn(12,2.5); // Returns: 10
494// p = quantdn([12,13,13.1,14,14.1,15,16],4); // Returns: [12,12,12,12,12,12,16]
495// q = quantdn([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,9,9,12]
496// r = quantdn([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[9,9,12]]
497function quantdn(x,y) =
498 assert( is_finite(y) && y>0, "The quantum `y` must be a positive value.")
499 is_num(x) ? floor(x/y)*y
500 : _floorall(x/y)*y;
501
502function _floorall(data) =
503 [for(x=data) is_list(x) ? _floorall(x) : floor(x)];
504
505
506// Function: quantup()
507// Synopsis: Returns `x` quantized uo to an integer multiple of `y`.
508// Topics: Math, Quantization
509// See Also: quant(), quantdn(), quantup()
510// Usage:
511// num = quantup(x, y);
512// Description:
513// Quantize a value `x` to an integer multiple of `y`, rounding up to the next multiple.
514// The value of `y` does NOT have to be an integer. If `x` is a list, then every item in
515// that list will be recursively quantized up.
516// Arguments:
517// x = The value or list to quantize.
518// y = Positive quantum to quantize to.
519// Example:
520// a = quantup(12,4); // Returns: 12
521// b = quantup(13,4); // Returns: 16
522// c = quantup(13.1,4); // Returns: 16
523// d = quantup(14,4); // Returns: 16
524// e = quantup(14.1,4); // Returns: 16
525// f = quantup(15,4); // Returns: 16
526// g = quantup(16,4); // Returns: 16
527// h = quantup(9,3); // Returns: 9
528// i = quantup(10,3); // Returns: 12
529// j = quantup(10.4,3); // Returns: 12
530// k = quantup(10.5,3); // Returns: 12
531// l = quantup(11,3); // Returns: 12
532// m = quantup(12,3); // Returns: 12
533// n = quantdn(11,2.5); // Returns: 12.5
534// o = quantdn(12,2.5); // Returns: 12.5
535// p = quantup([12,13,13.1,14,14.1,15,16],4); // Returns: [12,16,16,16,16,16,16]
536// q = quantup([9,10,10.4,10.5,11,12],3); // Returns: [9,12,12,12,12,12]
537// r = quantup([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,12,12],[12,12,12]]
538function quantup(x,y) =
539 assert( is_finite(y) && y>0, "The quantum `y` must be a positive value.")
540 is_num(x) ? ceil(x/y)*y
541 : _ceilall(x/y)*y;
542
543function _ceilall(data) =
544 [for(x=data) is_list(x) ? _ceilall(x) : ceil(x)];
545
546
547// Section: Constraints and Modulos
548
549// Function: constrain()
550// Synopsis: Returns a value constrained between `minval` and `maxval`, inclusive.
551// Topics: Math
552// See Also: posmod(), modang()
553// Usage:
554// val = constrain(v, minval, maxval);
555// Description:
556// Constrains value to a range of values between minval and maxval, inclusive.
557// Arguments:
558// v = value to constrain.
559// minval = minimum value to return, if out of range.
560// maxval = maximum value to return, if out of range.
561// Example:
562// a = constrain(-5, -1, 1); // Returns: -1
563// b = constrain(5, -1, 1); // Returns: 1
564// c = constrain(0.3, -1, 1); // Returns: 0.3
565// d = constrain(9.1, 0, 9); // Returns: 9
566// e = constrain(-0.1, 0, 9); // Returns: 0
567function constrain(v, minval, maxval) =
568 assert( is_finite(v+minval+maxval), "Input must be finite number(s).")
569 min(maxval, max(minval, v));
570
571
572// Function: posmod()
573// Synopsis: Returns the positive modulo of a value.
574// Topics: Math
575// See Also: constrain(), posmod(), modang()
576// Usage:
577// mod = posmod(x, m)
578// Description:
579// Returns the positive modulo `m` of `x`. Value returned will be in the range 0 ... `m`-1.
580// Arguments:
581// x = The value to constrain.
582// m = Modulo value.
583// Example:
584// a = posmod(-700,360); // Returns: 340
585// b = posmod(-270,360); // Returns: 90
586// c = posmod(-120,360); // Returns: 240
587// d = posmod(120,360); // Returns: 120
588// e = posmod(270,360); // Returns: 270
589// f = posmod(700,360); // Returns: 340
590// g = posmod(3,2.5); // Returns: 0.5
591function posmod(x,m) =
592 assert( is_finite(x) && is_finite(m) && !approx(m,0) , "Input must be finite numbers. The divisor cannot be zero.")
593 (x%m+m)%m;
594
595
596// Function: modang()
597// Synopsis: Returns an angle normalized to between -180º and 180º.
598// Topics: Math
599// See Also: constrain(), posmod(), modang()
600// Usage:
601// ang = modang(x);
602// Description:
603// Takes an angle in degrees and normalizes it to an equivalent angle value between -180 and 180.
604// Example:
605// a1 = modang(-700); // Returns: 20
606// a2 = modang(-270); // Returns: 90
607// a3 = modang(-120); // Returns: -120
608// a4 = modang(120); // Returns: 120
609// a5 = modang(270); // Returns: -90
610// a6 = modang(700); // Returns: -20
611function modang(x) =
612 assert( is_finite(x), "Input must be a finite number.")
613 let(xx = posmod(x,360)) xx<180? xx : xx-360;
614
615
616
617// Section: Operations on Lists (Sums, Mean, Products)
618
619// Function: sum()
620// Synopsis: Returns the sum of a list of values.
621// Topics: Math
622// See Also: mean(), median(), product(), cumsum()
623// Usage:
624// x = sum(v, [dflt]);
625// Description:
626// Returns the sum of all entries in the given consistent list.
627// If passed an array of vectors, returns the sum the vectors.
628// If passed an array of matrices, returns the sum of the matrices.
629// If passed an empty list, the value of `dflt` will be returned.
630// Arguments:
631// v = The list to get the sum of.
632// dflt = The default value to return if `v` is an empty list. Default: 0
633// Example:
634// sum([1,2,3]); // returns 6.
635// sum([[1,2,3], [3,4,5], [5,6,7]]); // returns [9, 12, 15]
636function sum(v, dflt=0) =
637 v==[]? dflt :
638 assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
639 is_finite(v[0]) || is_vector(v[0]) ? [for(i=v) 1]*v :
640 _sum(v,v[0]*0);
641
642function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
643
644
645
646
647// Function: mean()
648// Synopsis: Returns the mean value of a list of values.
649// Topics: Math, Statistics
650// See Also: sum(), mean(), median(), product()
651// Usage:
652// x = mean(v);
653// Description:
654// Returns the arithmetic mean/average of all entries in the given array.
655// If passed a list of vectors, returns a vector of the mean of each part.
656// Arguments:
657// v = The list of values to get the mean of.
658// Example:
659// mean([2,3,4]); // returns 3.
660// mean([[1,2,3], [3,4,5], [5,6,7]]); // returns [3, 4, 5]
661function mean(v) =
662 assert(is_list(v) && len(v)>0, "Invalid list.")
663 sum(v)/len(v);
664
665
666
667// Function: median()
668// Synopsis: Returns the median value of a list of values.
669// Topics: Math, Statistics
670// See Also: sum(), mean(), median(), product()
671// Usage:
672// middle = median(v)
673// Description:
674// Returns the median of the given vector.
675function median(v) =
676 assert(is_vector(v), "Input to median must be a vector")
677 len(v)%2 ? max( list_smallest(v, ceil(len(v)/2)) ) :
678 let( lowest = list_smallest(v, len(v)/2 + 1),
679 max = max(lowest),
680 imax = search(max,lowest,1),
681 max2 = max([for(i=idx(lowest)) if(i!=imax[0]) lowest[i] ])
682 )
683 (max+max2)/2;
684
685
686// Function: deltas()
687// Synopsis: Returns the deltas between a list of values.
688// Topics: Math, Statistics
689// See Also: sum(), mean(), median(), product()
690// Usage:
691// delts = deltas(v,[wrap]);
692// Description:
693// Returns a list with the deltas of adjacent entries in the given list, optionally wrapping back to the front.
694// The list should be a consistent list of numeric components (numbers, vectors, matrix, etc).
695// Given [a,b,c,d], returns [b-a,c-b,d-c].
696// Arguments:
697// v = The list to get the deltas of.
698// wrap = If true, wrap back to the start from the end. ie: return the difference between the last and first items as the last delta. Default: false
699// Example:
700// deltas([2,5,9,17]); // returns [3,4,8].
701// deltas([[1,2,3], [3,6,8], [4,8,11]]); // returns [[2,4,5], [1,2,3]]
702function deltas(v, wrap=false) =
703 assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.")
704 [for (p=pair(v,wrap)) p[1]-p[0]] ;
705
706
707// Function: cumsum()
708// Synopsis: Returns the running cumulative sum of a list of values.
709// Topics: Math, Statistics
710// See Also: sum(), mean(), median(), product()
711// Usage:
712// sums = cumsum(v);
713// Description:
714// Returns a list where each item is the cumulative sum of all items up to and including the corresponding entry in the input list.
715// If passed an array of vectors, returns a list of cumulative vectors sums.
716// Arguments:
717// v = The list to get the sum of.
718// Example:
719// cumsum([1,1,1]); // returns [1,2,3]
720// cumsum([2,2,2]); // returns [2,4,6]
721// cumsum([1,2,3]); // returns [1,3,6]
722// cumsum([[1,2,3], [3,4,5], [5,6,7]]); // returns [[1,2,3], [4,6,8], [9,12,15]]
723function cumsum(v) =
724 assert(is_consistent(v), "The input is not consistent." )
725 len(v)<=1 ? v :
726 _cumsum(v,_i=1,_acc=[v[0]]);
727
728function _cumsum(v,_i=0,_acc=[]) =
729 _i>=len(v) ? _acc :
730 _cumsum( v, _i+1, [ each _acc, _acc[len(_acc)-1] + v[_i] ] );
731
732
733
734// Function: product()
735// Synopsis: Returns the multiplicative product of a list of values.
736// Topics: Math, Statistics
737// See Also: sum(), mean(), median(), product(), cumsum()
738// Usage:
739// x = product(v);
740// Description:
741// Returns the product of all entries in the given list.
742// If passed a list of vectors of same dimension, returns a vector of products of each part.
743// If passed a list of square matrices, returns the resulting product matrix.
744// Arguments:
745// v = The list to get the product of.
746// Example:
747// product([2,3,4]); // returns 24.
748// product([[1,2,3], [3,4,5], [5,6,7]]); // returns [15, 48, 105]
749function product(v) =
750 assert( is_vector(v) || is_matrix(v) || ( is_matrix(v[0],square=true) && is_consistent(v)),
751 "Invalid input.")
752 _product(v, 1, v[0]);
753
754function _product(v, i=0, _tot) =
755 i>=len(v) ? _tot :
756 _product( v,
757 i+1,
758 ( is_vector(v[i])? v_mul(_tot,v[i]) : _tot*v[i] ) );
759
760
761
762// Function: cumprod()
763// Synopsis: Returns the running cumulative product of a list of values.
764// Topics: Math, Statistics
765// See Also: sum(), mean(), median(), product(), cumsum()
766// Usage:
767// prod_list = cumprod(list, [right]);
768// Description:
769// Returns a list where each item is the cumulative product of all items up to and including the corresponding entry in the input list.
770// If passed an array of vectors, returns a list of elementwise vector products. If passed a list of square matrices by default returns matrix
771// products multiplying on the left, so a list `[A,B,C]` will produce the output `[A,BA,CBA]`. If you set `right=true` then it returns
772// the product of multiplying on the right, so a list `[A,B,C]` will produce the output `[A,AB,ABC]` in that case.
773// Arguments:
774// list = The list to get the cumulative product of.
775// right = if true multiply matrices on the right
776// Example:
777// cumprod([1,3,5]); // returns [1,3,15]
778// cumprod([2,2,2]); // returns [2,4,8]
779// cumprod([[1,2,3], [3,4,5], [5,6,7]])); // returns [[1, 2, 3], [3, 8, 15], [15, 48, 105]]
780function cumprod(list,right=false) =
781 is_vector(list) ? _cumprod(list) :
782 assert(is_consistent(list), "Input must be a consistent list of scalars, vectors or square matrices")
783 assert(is_bool(right))
784 is_matrix(list[0]) ? assert(len(list[0])==len(list[0][0]), "Matrices must be square") _cumprod(list,right)
785 : _cumprod_vec(list);
786
787function _cumprod(v,right,_i=0,_acc=[]) =
788 _i==len(v) ? _acc :
789 _cumprod(
790 v, right, _i+1,
791 concat(
792 _acc,
793 [
794 _i==0 ? v[_i]
795 : right? _acc[len(_acc)-1]*v[_i]
796 : v[_i]*_acc[len(_acc)-1]
797 ]
798 )
799 );
800
801function _cumprod_vec(v,_i=0,_acc=[]) =
802 _i==len(v) ? _acc :
803 _cumprod_vec(
804 v, _i+1,
805 concat(
806 _acc,
807 [_i==0 ? v[_i] : v_mul(_acc[len(_acc)-1],v[_i])]
808 )
809 );
810
811
812
813// Function: convolve()
814// Synopsis: Returns the convolution of `p` and `q`.
815// Topics: Math, Statistics
816// See Also: sum(), mean(), median(), product(), cumsum()
817// Usage:
818// x = convolve(p,q);
819// Description:
820// Given two vectors, or one vector and a path or
821// two paths of the same dimension, finds the convolution of them.
822// If both parameter are vectors, returns the vector convolution.
823// If one parameter is a vector and the other a path,
824// convolves using products by scalars and returns a path.
825// If both parameters are paths, convolve using scalar products
826// and returns a vector.
827// The returned vector or path has length len(p)+len(q)-1.
828// Arguments:
829// p = The first vector or path.
830// q = The second vector or path.
831// Example:
832// a = convolve([1,1],[1,2,1]); // Returns: [1,3,3,1]
833// b = convolve([1,2,3],[1,2,1])); // Returns: [1,4,8,8,3]
834// c = convolve([[1,1],[2,2],[3,1]],[1,2,1])); // Returns: [[1,1],[4,4],[8,6],[8,4],[3,1]]
835// d = convolve([[1,1],[2,2],[3,1]],[[1,2],[2,1]])); // Returns: [3,9,11,7]
836function convolve(p,q) =
837 p==[] || q==[] ? [] :
838 assert( (is_vector(p) || is_matrix(p))
839 && ( is_vector(q) || (is_matrix(q) && ( !is_vector(p[0]) || (len(p[0])==len(q[0])) ) ) ) ,
840 "The inputs should be vectors or paths all of the same dimension.")
841 let( n = len(p),
842 m = len(q))
843 [for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) )
844 sum([for(j=[k1:k2]) p[i-j]*q[j] ])
845 ];
846
847
848
849// Function: sum_of_sines()
850// Synopsis: Returns the sum of one or more sine waves at a given angle.
851// Topics: Math, Statistics
852// See Also: sum(), mean(), median(), product(), cumsum()
853// Usage:
854// sum_of_sines(a,sines)
855// Description:
856// Given a list of sine waves, returns the sum of the sines at the given angle.
857// Each sine wave is given as an `[AMPLITUDE, FREQUENCY, PHASE_ANGLE]` triplet.
858// - `AMPLITUDE` is the height of the sine wave above (and below) `0`.
859// - `FREQUENCY` is the number of times the sine wave repeats in 360º.
860// - `PHASE_ANGLE` is the offset in degrees of the sine wave.
861// Arguments:
862// a = Angle to get the value for.
863// sines = List of [amplitude, frequency, phase_angle] items, where the frequency is the number of times the cycle repeats around the circle.
864// Example:
865// v = sum_of_sines(30, [[10,3,0], [5,5.5,60]]);
866function sum_of_sines(a, sines) =
867 assert( is_finite(a) && is_matrix(sines,undef,3), "Invalid input.")
868 sum([ for (s = sines)
869 let(
870 ss=point3d(s),
871 v=ss[0]*sin(a*ss[1]+ss[2])
872 ) v
873 ]);
874
875
876
877// Section: Random Number Generation
878
879// Function: rand_int()
880// Synopsis: Returns a random integer.
881// Topics: Random
882// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points(), exponential_rands()
883// Usage:
884// rand_int(minval, maxval, n, [seed]);
885// Description:
886// Return a list of random integers in the range of minval to maxval, inclusive.
887// Arguments:
888// minval = Minimum integer value to return.
889// maxval = Maximum integer value to return.
890// N = Number of random integers to return.
891// seed = If given, sets the random number seed.
892// Example:
893// ints = rand_int(0,100,3);
894// int = rand_int(-10,10,1)[0];
895function rand_int(minval, maxval, n, seed=undef) =
896 assert( is_finite(minval+maxval+n) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.")
897 assert(maxval >= minval, "Max value cannot be smaller than minval")
898 let (rvect = is_def(seed) ? rands(minval,maxval+1,n,seed) : rands(minval,maxval+1,n))
899 [for(entry = rvect) floor(entry)];
900
901
902// Function: random_points()
903// Synopsis: Returns a list of random points.
904// Topics: Random, Points
905// See Also: rand_int(), random_points(), random_polygon(), spherical_random_points()
906// Usage:
907// points = random_points(n, dim, [scale], [seed]);
908// Description:
909// Generate `n` uniform random points of dimension `dim` with data ranging from -scale to +scale.
910// The `scale` may be a number, in which case the random data lies in a cube,
911// or a vector with dimension `dim`, in which case each dimension has its own scale.
912// Arguments:
913// n = number of points to generate. Default: 1
914// dim = dimension of the points. Default: 2
915// scale = the scale of the point coordinates. Default: 1
916// seed = an optional seed for the random generation.
917function random_points(n, dim, scale=1, seed) =
918 assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.")
919 assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.")
920 assert( is_finite(scale) || is_vector(scale,dim), "The scale should be a number or a vector with length equal to d.")
921 let(
922 rnds = is_undef(seed)
923 ? rands(-1,1,n*dim)
924 : rands(-1,1,n*dim, seed) )
925 is_num(scale)
926 ? scale*[for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ]
927 : [for(i=[0:1:n-1]) [for(j=[0:dim-1]) scale[j]*rnds[i*dim+j] ] ];
928
929
930// Function: gaussian_rands()
931// Synopsis: Returns a list of random numbers with a gaussian distribution.
932// Topics: Random, Statistics
933// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points(), exponential_rands()
934// Usage:
935// arr = gaussian_rands([n],[mean], [cov], [seed]);
936// Description:
937// Returns a random number or vector with a Gaussian/normal distribution.
938// Arguments:
939// n = the number of points to return. Default: 1
940// mean = The average of the random value (a number or vector). Default: 0
941// cov = covariance matrix of the random numbers, or variance in the 1D case. Default: 1
942// seed = If given, sets the random number seed.
943function gaussian_rands(n=1, mean=0, cov=1, seed=undef) =
944 assert(is_num(mean) || is_vector(mean))
945 let(
946 dim = is_num(mean) ? 1 : len(mean)
947 )
948 assert((dim==1 && is_num(cov)) || is_matrix(cov,dim,dim),"mean and covariance matrix not compatible")
949 assert(is_undef(seed) || is_finite(seed))
950 let(
951 nums = is_undef(seed)? rands(0,1,dim*n*2) : rands(0,1,dim*n*2,seed),
952 rdata = [for (i = count(dim*n,0,2)) sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])]
953 )
954 dim==1 ? add_scalar(sqrt(cov)*rdata,mean) :
955 assert(is_matrix_symmetric(cov),"Supplied covariance matrix is not symmetric")
956 let(
957 L = cholesky(cov)
958 )
959 assert(is_def(L), "Supplied covariance matrix is not positive definite")
960 move(mean,list_to_matrix(rdata,dim)*transpose(L));
961
962
963// Function: exponential_rands()
964// Synopsis: Returns a list of random numbers with an exponential distribution.
965// Topics: Random, Statistics
966// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points()
967// Usage:
968// arr = exponential_rands([n], [lambda], [seed])
969// Description:
970// Returns random numbers with an exponential distribution with parameter lambda, and hence mean 1/lambda.
971// Arguments:
972// n = number of points to return. Default: 1
973// lambda = distribution parameter. The mean will be 1/lambda. Default: 1
974function exponential_rands(n=1, lambda=1, seed) =
975 assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
976 assert( is_num(lambda) && lambda>0, "The lambda parameter must be a positive number.")
977 let(
978 unif = is_def(seed) ? rands(0,1,n,seed=seed) : rands(0,1,n)
979 )
980 -(1/lambda) * [for(x=unif) x==1 ? 708.3964185322641 : ln(1-x)]; // Use ln(min_float) when x is 1
981
982
983// Function: spherical_random_points()
984// Synopsis: Returns a list of random points on the surface of a sphere.
985// Topics: Random, Points
986// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points()
987// Usage:
988// points = spherical_random_points([n], [radius], [seed]);
989// Description:
990// Generate `n` 3D uniformly distributed random points lying on a sphere centered at the origin with radius equal to `radius`.
991// Arguments:
992// n = number of points to generate. Default: 1
993// radius = the sphere radius. Default: 1
994// seed = an optional seed for the random generation.
995
996// See https://mathworld.wolfram.com/SpherePointPicking.html
997function spherical_random_points(n=1, radius=1, seed) =
998 assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
999 assert( is_num(radius) && radius>0, "The radius should be a non-negative number.")
1000 let( theta = is_undef(seed)
1001 ? rands(0,360,n)
1002 : rands(0,360,n, seed),
1003 cosphi = rands(-1,1,n))
1004 [for(i=[0:1:n-1]) let(
1005 sin_phi=sqrt(1-cosphi[i]*cosphi[i])
1006 )
1007 radius*[sin_phi*cos(theta[i]),sin_phi*sin(theta[i]), cosphi[i]]];
1008
1009
1010
1011// Function: random_polygon()
1012// Synopsis: Returns the CCW path of a simple random polygon.
1013// Topics: Random, Polygon
1014// See Also: random_points(), spherical_random_points()
1015// Usage:
1016// points = random_polygon([n], [size], [seed]);
1017// Description:
1018// Generate the `n` vertices of a random counter-clockwise simple 2d polygon
1019// inside a circle centered at the origin with radius `size`.
1020// Arguments:
1021// n = number of vertices of the polygon. Default: 3
1022// size = the radius of a circle centered at the origin containing the polygon. Default: 1
1023// seed = an optional seed for the random generation.
1024function random_polygon(n=3,size=1, seed) =
1025 assert( is_int(n) && n>2, "Improper number of polygon vertices.")
1026 assert( is_num(size) && size>0, "Improper size.")
1027 let(
1028 seed = is_undef(seed) ? rands(0,1,1)[0] : seed,
1029 cumm = cumsum(rands(0.1,10,n+1,seed)),
1030 angs = 360*cumm/cumm[n-1],
1031 rads = rands(.01,size,n,seed+cumm[0])
1032 )
1033 [for(i=count(n)) rads[i]*[cos(angs[i]), sin(angs[i])] ];
1034
1035
1036
1037// Section: Calculus
1038
1039// Function: deriv()
1040// Synopsis: Returns the first derivative estimate of a list of data.
1041// Topics: Math, Calculus
1042// See Also: deriv(), deriv2(), deriv3()
1043// Usage:
1044// x = deriv(data, [h], [closed])
1045// Description:
1046// Computes a numerical derivative estimate of the data, which may be scalar or vector valued.
1047// The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly.
1048// If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
1049// data[len(data)-1]. This function uses a symetric derivative approximation
1050// for internal points, f'(t) = (f(t+h)-f(t-h))/2h. For the endpoints (when closed=false) the algorithm
1051// uses a two point method if sufficient points are available: f'(t) = (3*(f(t+h)-f(t)) - (f(t+2*h)-f(t+h)))/2h.
1052// .
1053// If `h` is a vector then it is assumed to be nonuniform, with h[i] giving the sampling distance
1054// between data[i+1] and data[i], and the data values will be linearly resampled at each corner
1055// to produce a uniform spacing for the derivative estimate. At the endpoints a single point method
1056// is used: f'(t) = (f(t+h)-f(t))/h.
1057// Arguments:
1058// data = the list of the elements to compute the derivative of.
1059// h = the parametric sampling of the data.
1060// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
1061function deriv(data, h=1, closed=false) =
1062 assert( is_consistent(data) , "Input list is not consistent or not numerical.")
1063 assert( len(data)>=2, "Input `data` should have at least 2 elements.")
1064 assert( is_finite(h) || is_vector(h), "The sampling `h` must be a number or a list of numbers." )
1065 assert( is_num(h) || len(h) == len(data)-(closed?0:1),
1066 str("Vector valued `h` must have length ",len(data)-(closed?0:1)))
1067 is_vector(h) ? _deriv_nonuniform(data, h, closed=closed) :
1068 let( L = len(data) )
1069 closed
1070 ? [
1071 for(i=[0:1:L-1])
1072 (data[(i+1)%L]-data[(L+i-1)%L])/2/h
1073 ]
1074 : let(
1075 first = L<3 ? data[1]-data[0] :
1076 3*(data[1]-data[0]) - (data[2]-data[1]),
1077 last = L<3 ? data[L-1]-data[L-2]:
1078 (data[L-3]-data[L-2])-3*(data[L-2]-data[L-1])
1079 )
1080 [
1081 first/2/h,
1082 for(i=[1:1:L-2]) (data[i+1]-data[i-1])/2/h,
1083 last/2/h
1084 ];
1085
1086
1087function _dnu_calc(f1,fc,f2,h1,h2) =
1088 let(
1089 f1 = h2<h1 ? lerp(fc,f1,h2/h1) : f1 ,
1090 f2 = h1<h2 ? lerp(fc,f2,h1/h2) : f2
1091 )
1092 (f2-f1) / 2 / min(h1,h2);
1093
1094
1095function _deriv_nonuniform(data, h, closed) =
1096 let( L = len(data) )
1097 closed
1098 ? [for(i=[0:1:L-1])
1099 _dnu_calc(data[(L+i-1)%L], data[i], data[(i+1)%L], select(h,i-1), h[i]) ]
1100 : [
1101 (data[1]-data[0])/h[0],
1102 for(i=[1:1:L-2]) _dnu_calc(data[i-1],data[i],data[i+1], h[i-1],h[i]),
1103 (data[L-1]-data[L-2])/h[L-2]
1104 ];
1105
1106
1107// Function: deriv2()
1108// Synopsis: Returns the second derivative estimate of a list of data.
1109// Topics: Math, Calculus
1110// See Also: deriv(), deriv2(), deriv3()
1111// Usage:
1112// x = deriv2(data, [h], [closed])
1113// Description:
1114// Computes a numerical estimate of the second derivative of the data, which may be scalar or vector valued.
1115// The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly.
1116// If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
1117// data[len(data)-1]. For internal points this function uses the approximation
1118// f''(t) = (f(t-h)-2*f(t)+f(t+h))/h^2. For the endpoints (when closed=false),
1119// when sufficient points are available, the method is either the four point expression
1120// f''(t) = (2*f(t) - 5*f(t+h) + 4*f(t+2*h) - f(t+3*h))/h^2 or
1121// f''(t) = (35*f(t) - 104*f(t+h) + 114*f(t+2*h) - 56*f(t+3*h) + 11*f(t+4*h)) / 12h^2
1122// if five points are available.
1123// Arguments:
1124// data = the list of the elements to compute the derivative of.
1125// h = the constant parametric sampling of the data.
1126// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
1127function deriv2(data, h=1, closed=false) =
1128 assert( is_consistent(data) , "Input list is not consistent or not numerical.")
1129 assert( is_finite(h), "The sampling `h` must be a number." )
1130 let( L = len(data) )
1131 assert( L>=3, "Input list has less than 3 elements.")
1132 closed
1133 ? [
1134 for(i=[0:1:L-1])
1135 (data[(i+1)%L]-2*data[i]+data[(L+i-1)%L])/h/h
1136 ]
1137 :
1138 let(
1139 first =
1140 L==3? data[0] - 2*data[1] + data[2] :
1141 L==4? 2*data[0] - 5*data[1] + 4*data[2] - data[3] :
1142 (35*data[0] - 104*data[1] + 114*data[2] - 56*data[3] + 11*data[4])/12,
1143 last =
1144 L==3? data[L-1] - 2*data[L-2] + data[L-3] :
1145 L==4? -2*data[L-1] + 5*data[L-2] - 4*data[L-3] + data[L-4] :
1146 (35*data[L-1] - 104*data[L-2] + 114*data[L-3] - 56*data[L-4] + 11*data[L-5])/12
1147 ) [
1148 first/h/h,
1149 for(i=[1:1:L-2]) (data[i+1]-2*data[i]+data[i-1])/h/h,
1150 last/h/h
1151 ];
1152
1153
1154// Function: deriv3()
1155// Synopsis: Returns the third derivative estimate of a list of data.
1156// Topics: Math, Calculus
1157// See Also: deriv(), deriv2(), deriv3()
1158// Usage:
1159// x = deriv3(data, [h], [closed])
1160// Description:
1161// Computes a numerical third derivative estimate of the data, which may be scalar or vector valued.
1162// The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly.
1163// If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
1164// data[len(data)-1]. This function uses a five point derivative estimate, so the input data must include
1165// at least five points:
1166// f'''(t) = (-f(t-2*h)+2*f(t-h)-2*f(t+h)+f(t+2*h)) / 2h^3. At the first and second points from the end
1167// the estimates are f'''(t) = (-5*f(t)+18*f(t+h)-24*f(t+2*h)+14*f(t+3*h)-3*f(t+4*h)) / 2h^3 and
1168// f'''(t) = (-3*f(t-h)+10*f(t)-12*f(t+h)+6*f(t+2*h)-f(t+3*h)) / 2h^3.
1169// Arguments:
1170// data = the list of the elements to compute the derivative of.
1171// h = the constant parametric sampling of the data.
1172// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
1173function deriv3(data, h=1, closed=false) =
1174 assert( is_consistent(data) , "Input list is not consistent or not numerical.")
1175 assert( len(data)>=5, "Input list has less than 5 elements.")
1176 assert( is_finite(h), "The sampling `h` must be a number." )
1177 let(
1178 L = len(data),
1179 h3 = h*h*h
1180 )
1181 closed? [
1182 for(i=[0:1:L-1])
1183 (-data[(L+i-2)%L]+2*data[(L+i-1)%L]-2*data[(i+1)%L]+data[(i+2)%L])/2/h3
1184 ] :
1185 let(
1186 first=(-5*data[0]+18*data[1]-24*data[2]+14*data[3]-3*data[4])/2,
1187 second=(-3*data[0]+10*data[1]-12*data[2]+6*data[3]-data[4])/2,
1188 last=(5*data[L-1]-18*data[L-2]+24*data[L-3]-14*data[L-4]+3*data[L-5])/2,
1189 prelast=(3*data[L-1]-10*data[L-2]+12*data[L-3]-6*data[L-4]+data[L-5])/2
1190 ) [
1191 first/h3,
1192 second/h3,
1193 for(i=[2:1:L-3]) (-data[i-2]+2*data[i-1]-2*data[i+1]+data[i+2])/2/h3,
1194 prelast/h3,
1195 last/h3
1196 ];
1197
1198
1199// Section: Complex Numbers
1200
1201
1202// Function: complex()
1203// Synopsis: Replaces scalars in a list or matrix with complex number 2-vectors.
1204// Topics: Math, Complex Numbers
1205// See Also: c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1206// Usage:
1207// z = complex(list)
1208// Description:
1209// Converts a real valued number, vector or matrix into its complex analog
1210// by replacing all entries with a 2-vector that has zero imaginary part.
1211function complex(list) =
1212 is_num(list) ? [list,0] :
1213 [for(entry=list) is_num(entry) ? [entry,0] : complex(entry)];
1214
1215
1216// Function: c_mul()
1217// Synopsis: Multiplies two complex numbers.
1218// Topics: Math, Complex Numbers
1219// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1220// Usage:
1221// c = c_mul(z1,z2)
1222// Description:
1223// Multiplies two complex numbers, vectors or matrices, where complex numbers
1224// or entries are represented as vectors: [REAL, IMAGINARY]. Note that all
1225// entries in both arguments must be complex.
1226// Arguments:
1227// z1 = First complex number, vector or matrix
1228// z2 = Second complex number, vector or matrix
1229function c_mul(z1,z2) =
1230 is_matrix([z1,z2],2,2) ? _c_mul(z1,z2) :
1231 _combine_complex(_c_mul(_split_complex(z1), _split_complex(z2)));
1232
1233
1234function _split_complex(data) =
1235 is_vector(data,2) ? data
1236 : is_num(data[0][0]) ? [data*[1,0], data*[0,1]]
1237 : [
1238 [for(vec=data) vec * [1,0]],
1239 [for(vec=data) vec * [0,1]]
1240 ];
1241
1242
1243function _combine_complex(data) =
1244 is_vector(data,2) ? data
1245 : is_num(data[0][0]) ? [for(i=[0:len(data[0])-1]) [data[0][i],data[1][i]]]
1246 : [for(i=[0:1:len(data[0])-1])
1247 [for(j=[0:1:len(data[0][0])-1])
1248 [data[0][i][j], data[1][i][j]]]];
1249
1250
1251function _c_mul(z1,z2) =
1252 [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ];
1253
1254
1255// Function: c_div()
1256// Synopsis: Divides two complex numbers.
1257// Topics: Math, Complex Numbers
1258// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1259// Usage:
1260// x = c_div(z1,z2)
1261// Description:
1262// Divides two complex numbers represented by 2D vectors.
1263// Returns a complex number as a 2D vector [REAL, IMAGINARY].
1264// Arguments:
1265// z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
1266// z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
1267function c_div(z1,z2) =
1268 assert( is_vector(z1,2) && is_vector(z2), "Complex numbers should be represented by 2D vectors." )
1269 assert( !approx(z2,0), "The divisor `z2` cannot be zero." )
1270 let(den = z2.x*z2.x + z2.y*z2.y)
1271 [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x - z1.x*z2.y)/den];
1272
1273
1274// Function: c_conj()
1275// Synopsis: Returns the complex conjugate of the input.
1276// Topics: Math, Complex Numbers
1277// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1278// Usage:
1279// w = c_conj(z)
1280// Description:
1281// Computes the complex conjugate of the input, which can be a complex number,
1282// complex vector or complex matrix.
1283function c_conj(z) =
1284 is_vector(z,2) ? [z.x,-z.y] :
1285 [for(entry=z) c_conj(entry)];
1286
1287
1288// Function: c_real()
1289// Synopsis: Returns the real part of a complex number, vector or matrix..
1290// Topics: Math, Complex Numbers
1291// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1292// Usage:
1293// x = c_real(z)
1294// Description:
1295// Returns real part of a complex number, vector or matrix.
1296function c_real(z) =
1297 is_vector(z,2) ? z.x
1298 : is_num(z[0][0]) ? z*[1,0]
1299 : [for(vec=z) vec * [1,0]];
1300
1301
1302// Function: c_imag()
1303// Synopsis: Returns the imaginary part of a complex number, vector or matrix..
1304// Topics: Math, Complex Numbers
1305// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1306// Usage:
1307// x = c_imag(z)
1308// Description:
1309// Returns imaginary part of a complex number, vector or matrix.
1310function c_imag(z) =
1311 is_vector(z,2) ? z.y
1312 : is_num(z[0][0]) ? z*[0,1]
1313 : [for(vec=z) vec * [0,1]];
1314
1315
1316// Function: c_ident()
1317// Synopsis: Returns an n by n complex identity matrix.
1318// Topics: Math, Complex Numbers
1319// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1320// Usage:
1321// I = c_ident(n)
1322// Description:
1323// Produce an n by n complex identity matrix
1324function c_ident(n) = [for (i = [0:1:n-1]) [for (j = [0:1:n-1]) (i==j)?[1,0]:[0,0]]];
1325
1326
1327// Function: c_norm()
1328// Synopsis: Returns the norm of a complex number or vector.
1329// Topics: Math, Complex Numbers
1330// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
1331// Usage:
1332// n = c_norm(z)
1333// Description:
1334// Compute the norm of a complex number or vector.
1335function c_norm(z) = norm_fro(z);
1336
1337
1338// Section: Polynomials
1339
1340// Function: quadratic_roots()
1341// Synopsis: Computes roots for the quadratic equation.
1342// Topics: Math, Geometry, Complex Numbers
1343// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add()
1344// Usage:
1345// roots = quadratic_roots(a, b, c, [real])
1346// Description:
1347// Computes roots of the quadratic equation a*x^2+b*x+c==0, where the
1348// coefficients are real numbers. If real is true then returns only the
1349// real roots. Otherwise returns a pair of complex values. This method
1350// may be more reliable than the general root finder at distinguishing
1351// real roots from complex roots.
1352// Algorithm from: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
1353function quadratic_roots(a,b,c,real=false) =
1354 real ? [for(root = quadratic_roots(a,b,c,real=false)) if (root.y==0) root.x]
1355 :
1356 is_undef(b) && is_undef(c) && is_vector(a,3) ? quadratic_roots(a[0],a[1],a[2]) :
1357 assert(is_num(a) && is_num(b) && is_num(c))
1358 assert(a!=0 || b!=0 || c!=0, "Quadratic must have a nonzero coefficient")
1359 a==0 && b==0 ? [] : // No solutions
1360 a==0 ? [[-c/b,0]] :
1361 let(
1362 descrim = b*b-4*a*c,
1363 sqrt_des = sqrt(abs(descrim))
1364 )
1365 descrim < 0 ? // Complex case
1366 [[-b, sqrt_des],
1367 [-b, -sqrt_des]]/2/a :
1368 b<0 ? // b positive
1369 [[2*c/(-b+sqrt_des),0],
1370 [(-b+sqrt_des)/a/2,0]]
1371 : // b negative
1372 [[(-b-sqrt_des)/2/a, 0],
1373 [2*c/(-b-sqrt_des),0]];
1374
1375
1376// Function: polynomial()
1377// Synopsis: Evaluate a polynomial at a real or complex value.
1378// Topics: Math, Complex Numbers
1379// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
1380// Usage:
1381// x = polynomial(p, z)
1382// Description:
1383// Evaluates specified real polynomial, p, at the complex or real input value, z.
1384// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
1385// where a_n is the z^n coefficient. Polynomial coefficients are real.
1386// The result is a number if `z` is a number and a complex number otherwise.
1387function polynomial(p,z,k,total) =
1388 is_undef(k)
1389 ? assert( is_vector(p) , "Input polynomial coefficients must be a vector." )
1390 assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
1391 polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0])
1392 : k==len(p) ? total
1393 : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : c_mul(total,z)+[p[k],0]);
1394
1395
1396// Function: poly_mult()
1397// Synopsis: Compute product of two polynomials, returning a polynomial.
1398// Topics: Math
1399// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
1400// Usage:
1401// x = polymult(p,q)
1402// x = polymult([p1,p2,p3,...])
1403// Description:
1404// Given a list of polynomials represented as real algebraic coefficient lists, with the highest degree coefficient first,
1405// computes the coefficient list of the product polynomial.
1406function poly_mult(p,q) =
1407 is_undef(q) ?
1408 len(p)==2
1409 ? poly_mult(p[0],p[1])
1410 : poly_mult(p[0], poly_mult(list_tail(p)))
1411 :
1412 assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
1413 p*p==0 || q*q==0
1414 ? [0]
1415 : _poly_trim(convolve(p,q));
1416
1417
1418// Function: poly_div()
1419// Synopsis: Returns the polynomial quotient and remainder results of dividing two polynomials.
1420// Topics: Math
1421// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
1422// Usage:
1423// [quotient,remainder] = poly_div(n,d)
1424// Description:
1425// Computes division of the numerator polynomial by the denominator polynomial and returns
1426// a list of two polynomials, [quotient, remainder]. If the division has no remainder then
1427// the zero polynomial [0] is returned for the remainder. Similarly if the quotient is zero
1428// the returned quotient will be [0].
1429function poly_div(n,d) =
1430 assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
1431 let( d = _poly_trim(d),
1432 n = _poly_trim(n) )
1433 assert( d!=[0] , "Denominator cannot be a zero polynomial." )
1434 n==[0]
1435 ? [[0],[0]]
1436 : _poly_div(n,d,q=[]);
1437
1438function _poly_div(n,d,q) =
1439 len(n)<len(d) ? [q,_poly_trim(n)] :
1440 let(
1441 t = n[0] / d[0],
1442 newq = concat(q,[t]),
1443 newn = [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
1444 )
1445 _poly_div(newn,d,newq);
1446
1447
1448/// Internal Function: _poly_trim()
1449/// Usage:
1450/// _poly_trim(p, [eps])
1451/// Description:
1452/// Removes leading zero terms of a polynomial. By default zeros must be exact,
1453/// or give epsilon for approximate zeros. Returns [0] for a zero polynomial.
1454function _poly_trim(p,eps=0) =
1455 let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i])
1456 len(nz)==0 ? [0] : list_tail(p,nz[0]);
1457
1458
1459// Function: poly_add()
1460// Synopsis: Returns the polynomial sum of adding two polynomials.
1461// Topics: Math
1462// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
1463// Usage:
1464// sum = poly_add(p,q)
1465// Description:
1466// Computes the sum of two polynomials.
1467function poly_add(p,q) =
1468 assert( is_vector(p) && is_vector(q), "Invalid input polynomial(s)." )
1469 let( plen = len(p),
1470 qlen = len(q),
1471 long = plen>qlen ? p : q,
1472 short = plen>qlen ? q : p
1473 )
1474 _poly_trim(long + concat(repeat(0,len(long)-len(short)),short));
1475
1476
1477// Function: poly_roots()
1478// Synopsis: Returns all complex valued roots of the given real polynomial.
1479// Topics: Math, Complex Numbers
1480// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
1481// Usage:
1482// roots = poly_roots(p, [tol]);
1483// Description:
1484// Returns all complex roots of the specified real polynomial p.
1485// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
1486// where a_n is the z^n coefficient. The tol parameter gives
1487// the stopping tolerance for the iteration. The polynomial
1488// must have at least one non-zero coefficient. Convergence is poor
1489// if the polynomial has any repeated roots other than zero.
1490// Arguments:
1491// p = polynomial coefficients with higest power coefficient first
1492// tol = tolerance for iteration. Default: 1e-14
1493
1494// Uses the Aberth method https://en.wikipedia.org/wiki/Aberth_method
1495//
1496// Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996.
1497// https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method
1498function poly_roots(p,tol=1e-14,error_bound=false) =
1499 assert( is_vector(p), "Invalid polynomial." )
1500 let( p = _poly_trim(p,eps=0) )
1501 assert( p!=[0], "Input polynomial cannot be zero." )
1502 p[len(p)-1] == 0 ? // Strip trailing zero coefficients
1503 let( solutions = poly_roots(list_head(p),tol=tol, error_bound=error_bound))
1504 (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
1505 : [[0,0], each solutions]) :
1506 len(p)==1 ? (error_bound ? [[],[]] : []) : // Nonzero constant case has no solutions
1507 len(p)==2 ? let( solution = [[-p[1]/p[0],0]]) // Linear case needs special handling
1508 (error_bound ? [solution,[0]] : solution)
1509 :
1510 let(
1511 n = len(p)-1, // polynomial degree
1512 pderiv = [for(i=[0:n-1]) p[i]*(n-i)],
1513
1514 s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)], // Error bound polynomial from Bini
1515
1516 // Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf
1517 beta = -p[1]/p[0]/n,
1518 r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n),
1519 init = [for(i=[0:1:n-1]) // Initial guess for roots
1520 let(angle = 360*i/n+270/n/PI)
1521 [beta,0]+r*[cos(angle),sin(angle)]
1522 ],
1523 roots = _poly_roots(p,pderiv,s,init,tol=tol),
1524 error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) /
1525 abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0
1526 )
1527 error_bound ? [roots, error] : roots;
1528
1529// Internal function
1530// p = polynomial
1531// pderiv = derivative polynomial of p
1532// z = current guess for the roots
1533// tol = root tolerance
1534// i=iteration counter
1535function _poly_roots(p, pderiv, s, z, tol, i=0) =
1536 assert(i<45, str("Polyroot exceeded iteration limit. Current solution:", z))
1537 let(
1538 n = len(z),
1539 svals = [for(zk=z) tol*polynomial(s,norm(zk))],
1540 p_of_z = [for(zk=z) polynomial(p,zk)],
1541 done = [for(k=[0:n-1]) norm(p_of_z[k])<=svals[k]],
1542 newton = [for(k=[0:n-1]) c_div(p_of_z[k], polynomial(pderiv,z[k]))],
1543 zdiff = [for(k=[0:n-1]) sum([for(j=[0:n-1]) if (j!=k) c_div([1,0], z[k]-z[j])])],
1544 w = [for(k=[0:n-1]) done[k] ? [0,0] : c_div( newton[k],
1545 [1,0] - c_mul(newton[k], zdiff[k]))]
1546 )
1547 all(done) ? z : _poly_roots(p,pderiv,s,z-w,tol,i+1);
1548
1549
1550// Function: real_roots()
1551// Synopsis: Returns all real roots of the given real polynomial.
1552// Topics: Math, Complex Numbers
1553// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
1554// Usage:
1555// roots = real_roots(p, [eps], [tol])
1556// Description:
1557// Returns the real roots of the specified real polynomial p.
1558// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
1559// where a_n is the x^n coefficient. This function works by
1560// computing the complex roots and returning those roots where
1561// the imaginary part is closed to zero. By default it uses a computed
1562// error bound from the polynomial solver to decide whether imaginary
1563// parts are zero. You can specify eps, in which case the test is
1564// z.y/(1+norm(z)) < eps. Because
1565// of poor convergence and higher error for repeated roots, such roots may
1566// be missed by the algorithm because error can make their imaginary parts
1567// large enough to appear non-zero.
1568// Arguments:
1569// p = polynomial to solve as coefficient list, highest power term first
1570// eps = used to determine whether imaginary parts of roots are zero
1571// tol = tolerance for the complex polynomial root finder
1572
1573function real_roots(p,eps=undef,tol=1e-14) =
1574 assert( is_vector(p), "Invalid polynomial." )
1575 let( p = _poly_trim(p,eps=0) )
1576 assert( p!=[0], "Input polynomial cannot be zero." )
1577 let(
1578 roots_err = poly_roots(p,error_bound=true),
1579 roots = roots_err[0],
1580 err = roots_err[1]
1581 )
1582 is_def(eps)
1583 ? [for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x]
1584 : [for(i=idx(roots)) if (abs(roots[i].y)<=err[i]) roots[i].x];
1585
1586
1587// Section: Operations on Functions
1588
1589// Function: root_find()
1590// Synopsis: Finds a root of the given continuous function.
1591// Topics: Math
1592// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
1593// Usage:
1594// x = root_find(f, x0, x1, [tol])
1595// Description:
1596// Find a root of the continuous function f where the sign of f(x0) is different
1597// from the sign of f(x1). The function f is a function literal accepting one
1598// argument. You must have a version of OpenSCAD that supports function literals
1599// (2021.01 or newer). The tolerance (tol) specifies the accuracy of the solution:
1600// abs(f(x)) < tol * yrange, where yrange is the range of observed function values.
1601// This function can only find roots that cross the x axis: it cannot find the
1602// the root of x^2.
1603// Arguments:
1604// f = function literal for a scalar-valued single variable function
1605// x0 = endpoint of interval to search for root
1606// x1 = second endpoint of interval to search for root
1607// tol = tolerance for solution. Default: 1e-15
1608
1609// The algorithm is based on Brent's method and is a combination of
1610// bisection and inverse quadratic approximation, where bisection occurs
1611// at every step, with refinement using inverse quadratic approximation
1612// only when that approximation gives a good result. The detail
1613// of how to decide when to use the quadratic came from an article
1614// by Crenshaw on "The World's Best Root Finder".
1615// https://www.embedded.com/worlds-best-root-finder/
1616function root_find(f,x0,x1,tol=1e-15) =
1617 let(
1618 y0 = f(x0),
1619 y1 = f(x1),
1620 yrange = y0<y1 ? [y0,y1] : [y1,y0]
1621 )
1622 // Check endpoints
1623 y0==0 || _rfcheck(x0, y0,yrange,tol) ? x0 :
1624 y1==0 || _rfcheck(x1, y1,yrange,tol) ? x1 :
1625 assert(y0*y1<0, "Sign of function must be different at the interval endpoints")
1626 _rootfind(f,[x0,x1],[y0,y1],yrange,tol);
1627
1628function _rfcheck(x,y,range,tol) =
1629 assert(is_finite(y), str("Function not finite at ",x))
1630 abs(y) < tol*(range[1]-range[0]);
1631
1632// xpts and ypts are arrays whose first two entries contain the
1633// interval bracketing the root. Extra entries are ignored.
1634// yrange is the total observed range of y values (used for the
1635// tolerance test).
1636function _rootfind(f, xpts, ypts, yrange, tol, i=0) =
1637 assert(i<100, "root_find did not converge to a solution")
1638 let(
1639 xmid = (xpts[0]+xpts[1])/2,
1640 ymid = f(xmid),
1641 yrange = [min(ymid, yrange[0]), max(ymid, yrange[1])]
1642 )
1643 _rfcheck(xmid, ymid, yrange, tol) ? xmid :
1644 let(
1645 // Force root to be between x0 and midpoint
1646 y = ymid * ypts[0] < 0 ? [ypts[0], ymid, ypts[1]]
1647 : [ypts[1], ymid, ypts[0]],
1648 x = ymid * ypts[0] < 0 ? [xpts[0], xmid, xpts[1]]
1649 : [xpts[1], xmid, xpts[0]],
1650 v = y[2]*(y[2]-y[0]) - 2*y[1]*(y[1]-y[0])
1651 )
1652 v <= 0 ? _rootfind(f,x,y,yrange,tol,i+1) // Root is between first two points, extra 3rd point doesn't hurt
1653 :
1654 let( // Do quadratic approximation
1655 B = (x[1]-x[0]) / (y[1]-y[0]),
1656 C = y*[-1,2,-1] / (y[2]-y[1]) / (y[2]-y[0]),
1657 newx = x[0] - B * y[0] *(1-C*y[1]),
1658 newy = f(newx),
1659 new_yrange = [min(yrange[0],newy), max(yrange[1], newy)],
1660 // select interval that contains the root by checking sign
1661 yinterval = newy*y[0] < 0 ? [y[0],newy] : [newy,y[1]],
1662 xinterval = newy*y[0] < 0 ? [x[0],newx] : [newx,x[1]]
1663 )
1664 _rfcheck(newx, newy, new_yrange, tol)
1665 ? newx
1666 : _rootfind(f, xinterval, yinterval, new_yrange, tol, i+1);
1667
1668
1669
1670// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap